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A method for optimizing a central receiver solar thermal electric power plant is studied. We parametrize the plant design as a 
function of eleven design variables and reduce the problem of finding optimal designs to the numerical problem of finding the 
minimum of a function of several variables. This minimization problem is attacked with different algorithms both local and global 
in nature. We find that all algorithms find the same minimum of the objective function. The performance of each of the algorithms 
and the resulting designs are studied for two typical cases. 

We describe a method to evaluate the impact of design variables in the plant performance. This method will tell us what variables 
are key to the optimal plant design and which ones are less important. This information can be used to further improve the plant 
design and to accelerate the optimization procedure. 

Keywords: optimization; solar thermal electric plant design; field layout; collector field design. 



1. Introduction 

One of the main tasks in the conversion of solar energy into 
electricity by solar power plants is to work out an optimized 
plant design. In this type of plant, the energy collector subsys- 
tems (heliostats, field receivers) represent a very important part 
of the cost break-down structure. Therefore, the use of detailed 
computer programs is of great interest in order to optimize the 
plant design. 

The two main conceptual ingredients for a solar plant opti- 
mization code are: 
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1 . Reduction of the plant design to the value of certain design 
variables. 

2. An optimization criteria: This means having a function 
that computes the objective quantity (i.e. total annual 
power output, cost per produced power, etc. . . ) as a func- 
tion of the design variables. In general we will use the cost 
of the energy produced by the plant as the optimization cri- 
teria, and therefore use the terms optimize and minimize as 
interchangeable. 

After the plant design is completed it becomes crucial to 
understand what role each of the design variables play in the 
optimal design. Some variables can be slightly moved away 
from the optimal value without impacting the plant perfor- 
mance while others can not be changed without a severe im- 
pact in the plant performance. We will give a precise defini- 
tion of a quantity (we call it uncertainty) that will measure the 
importance of each variable in a plant design. We will give a 
precise mathematical definition of this uncertainty associated 
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with a variable, and show how to compute these quantities for 
a general plant design. 

The paper is organized as follows. In section [2] we will ex- 
plain our choice of plant design variables. Sections [3] and ^de- 
scribe additional information needed to perform the optimiza- 
tion. We perform the numerical optimization using different al- 
gorithms, both global and local. Section[5]describes in detail the 
optimization procedure and our choice of three different algo- 
rithms: first we use a fast local optimizer especially designed to 
solve this problem, second we use complex optimization library 
that include a Monte-Carlo search, with the potential ability to 
jump over function barriers and find global optima. Finally we 
use a genetic algorithm, generally used for difficult optimiza- 
tions and problems with a complex fitness landscape. 

Section [6] attack the very important problem of determining 
what role each of the design variables have in the plant perfor- 
mance. As we will see, the Hessian of the objective function 
evaluated at the minimum will provide information on how flat 
each of the directions in the minimum are. This information 
will not only give us the uncertainty associated with each design 
variable, but also help us deciding if the function have several 
local optima and if the three algorithms have found the same 
design. 

Finally, in section [7] we apply these optimization algorithms 
and analysis tools to two typical plant designs, and in section [8] 
we summarize our main results and present some perspectives 
for further studies. 



2. Description of plant design variables 

In choosing the variables that determine the layout of our 
plant, one has to take into account that an optimization pro- 
cedure requires multiple evaluations of the objective function. 
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Figure 1 : An example of the first two lines of heliostats in a field. The azimuthal 
angle starts in the north with the value 8 = 0. In green we have the heliostats 
that belong to the first line. The distance of this line to the tower is given by i?i . 
The distance of all the lines R„ is determined trough the recursion relation of 
Eq. [T] We see that this distance dependes on the azimuthal coordinate through 
the variable dg (Eq. [2}- The azimuthal distance between heliostats of the first 
line (i)J,D},...) are determined through the recursive relation Eq. |3] The az- 
imuthal position of the second and subsequent lines of heliostats is determined 
by positioning them at the average azimuthal position of the heliostats of the 
previous line (radial staggered) until a transition line is placed. 

Thus a proper choice of design variables must have the CPU 
cost of the optimization in mind. 

In this section we will present our choice for plant design 
variables. This choice has been made with the following things 
in mind 

• A plant should be circular-like to minimise the blocking 
and shadowing effect. 

• Consecutive files of heliostats should be allocated in radial 
staggered position to minimise blocking and shadowing 
effect, and keep a compact field. 

2.1. Collector field variables 

The distance of the n' h heliostat line to the tower (R„) is given 
by the linear recursion formulae 

R n = a + (fii + l)R n -i (1) 
that defines the two variables 

V01 ao : Initial spacing between heliostat rows [m], 
V02 a\ : Increasing spacing between heliostat rows. 

The recursion relation of Eq. Q] is solved by the condition 
= Rbase that specifies the distance between the tower and the 
first heliostat line. This space is usually used for operational 
purposes, like the administrative building, roads, etc. . . There 
is a minimal distance between heliostats R mm related with the 
size of the heliostats and some practical needs. R„ is set to the 
maximum between the value given by the recursion relation of 
Eq.Q]and R min . 

The shape of the lines of heliostats of an optimal field does 
not have to be perfectly circular. We consider the possibility 
of non circular shapes by adding to the previous radial distance 



an increment that dependes on the azimuthal position of the 
heliostat (we will call it 9). 

KR = / de ° for < 6» < 7T 

\ d e {2n - 8) for n < 8 < 2n ' K ' 

This defines another collector variable 

V03 dg : Correction of the radial distance with the azimuthal 
position. 

Intuitively when dg > the lines of heliostats will be closer 
to the tower in the north part of the field, whereas for dg < 
heliostats will be closer to the tower in the south. Usually 
optimal designs are more compact in the south (and hence dg < 
0) where there is less blocking and shadowing effect between 
heliostats. 

A proper optimal layout should not only find the optimum 
radial spacing between heliostats of the same line, but also de- 
termine the azimuthal distance between them. We will use a 
similar technique than the one used for the radial distance. We 
need to determine how this azimuthal distance depends on both 
the azimuthal angle and the radial distance. If we call the 
azimuthal distance for the heliostat that is just in the north (zero 
azimuth) in the first line of heliostats, and number the heliostats 
in the same line with the index a = 0, 1, . . ., the azimuthal dis- 
tance as a function of the azimuthal angle 8 is given by 

D l a = D l a _, + e e 8 (3) 

that, defines an additional variable 

V04 eg : Variation of the azimuthal distance with the azimuthal 
angle. 

To start this recursion relation one should give Dq as input. We 
will comment about this later. 

Variable ([V04]) determine the azimuthal position only of the 
first line of heliostats. The remaining lines of heliostats are sit- 
uated at the radial distance determined by the variables ([V01- 
V03]) and with radial staggered positions. This means that their 
azimuthal angle is the average of the azimuthal position of the 
two heliostats in front of it (see Fig . [TJ . 

This rule of formation will increase azimuthal spacing be- 
tween heliostats very fast. So it is convenient, after a certain 
amount of lines, to give an extra space and restart again the rule 
of formation forgetting the previous heliostats. This extra space 
increases with the distance to the tower. These re-starting lines 
are called transition lines, and heliostats between two transition 
lines are said to belong to the same group. The first line of he- 
liostats of each group is at a distance from the previous heliostat 
line given by 

A = (1 +ao + aiR)6 + e (4) 

where R is the last distance between lines given by Eq.Q] Eq.|4] 
defines two more optimization variables 

V05 e : Extra distance for transition lines. 
V06 5 : Increment of distance in transition lines with distance 
from the tower. 
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The variable e determines the extra spacing in a transition line 
to avoid large blocking and shadowing effects. Variable 5 sets 
an extra space that is proportional to the radial distance of the 
transition line. In Fig.|2]we can see the transition lines of a 3000 
heliostats plant. 
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Figure 2: Close look to the first transition lines of the field layout of an optimal 
3000 heliostat plant. We can see the transition lines, the number of heliostats in 
each these transition lines are 2-3-3-4-4-5-... (see text for more details). This 
optimal layout correspond to the output of the NSPOC code (see Sections [5] 
andUJfor more details). 

After a transition line we have to determine the azimuthal 
distance again of the next group of heliostats. This is done with 
the same recursion relation Eq. [3] 

K = DU+e e 6 (5) 

where the index a = 1,2... label groups of heliostats. As was 
commented before these recursion relations need D a Q as an ini- 
tial condition. These quantities are determined through 

Dg = (Z> + W (6) 

that essentially allow the azimuthal distance between heliostats 
in the north to increase/decrease with the radial distance of the 
group. This last recursion relation needs D l as an input. This 
is another optimization parameter completing the 8 parameters 
needed to optimize a field layout 

V07 b : Azimuthal distance dependence with the radial distance 
of the group. 

V08 : Initial azimuthal distance in the first line of heliostats 
of the plant. 

It is understood that there is a minimal azimuthal distance D m ,„ 
between heliostats that is given by the size of the heliostat plus 
some arbitrary distance needed for operational purposes. 

The only additional information needed is the number of he- 
liostat lines of each group. These are determined by trial and 
error. Reasonable results are obtained with an increasing se- 
quence like 2 -3 -3 -4-4-5 -5 -6 -6 -7 -7... (see Fig.|2] 
for an example). The conclusions of this paper are unchanged 
by the details like how many heliostats each group has, one 
only needs to fine tune this numbers when looking for the final 
design of a plant, and this can be easily done in an automated 
way. 



It is important to remark here that the set of variables ([V01- 
V08]) are compatible with much more conventional field de- 
signs. For example having no transition lines is easily achieved 
by adjusting e and 5. A constant azimuthal distance between 
heliostats by setting b — eg = 0, and a constant separation be- 
tween lines corresponds to the choice a\ — 0. Finally dg = 
give perfectly circular heliostat lines. 

The existence of some of the variables allow more complex 
designs that we think can improve the plant performance, but 
this complex design is not imposed. The optimization process 
will choose between the different options. 

Following these rules one generates a general field layout. 
We chose to fix the number of heliostats of our plant design 
(Nhei)- In order to choose these heliostats we pick up the ones 
which contribute with more power to the receiver. 

800 I 1 , 1 1 1 1 
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Figure 3: Field layout of an optimal 900 heliostat plant. This optimal layout 
correspond to the output of the NSPOC code (see Sections [5] and for more 
details). During the optimization more heliostats are generated but not selected 
for the final layout due to their worse performance compared to the selected 
ones (see Fig. [4). 
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Figure 4: Heliostats not selected for the final layout due to their worse perfor- 
mance. This image shows a subset of all the generated heliostats during the 
optimization (see text and Fig.fJ]for more details.) 

In Fig. [3] and Fig. [4] we can see how this procedure works 
in practice for a 900 heliostat plant. During the optimization 
process, we generate a field layout with far more than 900 he- 
liostats. To draw the final layout we simply pick the "best" he- 
liostats of the field (those shown in Fig. [3), but more heliostats 
are generated and rejected due to their worse performance (see 
Fig.©. 
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Figure 5: Field layout of a plant with a tower with three receivers. One of 
them points to the north, and the others at 120°. As can be seen the presence 
of two additional receivers increases dramatically the interception efficiency of 
heliostats at the east and west, making them to contribute with more power to 
the receiver than other heliostats in the north but more distant. This effect make 
the optimal field layout to have a clover shape. 



Note that this rule to select the heliostats does not necessar- 
ily impose a circular-like field layout. Simply this layout turns 
out to be optimal for the problems we are studying here. As 
an example of a non-circular optimal field, we show in Fig. 
the field layout of a more compl icated plant design with t hree 
cavity receivers in the tower (see lCrespo and Ramosl (120091) for 
more details). As can be seen in Fig. [5] having three cavity 
receivers dramatically improves the interception efficiency of 
some heliostats in the east and west part of the field, making 
this heliostats to perform better that heliostats that are far in 
the north. In this case the optimal field layout does not have a 
circular shape, but clover-like. 

2.2. Receiver variables 

We can design either a north field plant or a circular plant. In 
the first case we will consider that the receiver consists on an 
circular aperture^ pointing to the north that lets the reflected so- 
lar rays to enter. This receiver is characterized by the following 
variables 

V09 hj : Tower height [m], 

V10 r : Aperture radius [rad]. 

VI 1 ei : Aperture inclination [rad]. 

In the second case the receiver consists in a cylinder that can 
absorb radiation coming from all directions. In this case the 
receiver is characterized by its position and size0, parametrized 
by the following variables 

V09 h, : Tower height [to]. 
VI r : Receiver radius [rad]. 
VI 1 h r : Receiver height [to]. 



"We have considered more complex aperture shapes, like elliptical or rect- 
angular. Although some of this designs perform better than a circular aperture, 
the conclusions of the paper remains completely unaltered, and keep the dis- 
cussion simpler. 

3 In real world designs one should worry about the maximum power den- 
sity absorbed by the receiver, since this is constrained by material properties. 
This can (and should) be included in the optimization process, but we will not 
address this problem here. 



3. Additional input: minimization parameters 

To compute the full output of a solar power plant, we need 
some additional input related with the heliostat characteristics 
and plant location. These parameters are not treated as variables 
in the minimization process. 

These extra input are seven parameters that define the he- 
liostats characteristics and plant location, plus insolation data 
from the plant location. 

The study on how the plant performance depends on these 
parameters (what we call parametric analysis) is definitively a 
very interesting subject that can answer questions like What is 
the proper heliostat size? How much influence the optical qual- 
ity of an heliostat the plant performance?, etc... Nevertheless 
these questions are beyond the scope of the pr esent work and 
need further investigations (ICrespo et al.i 1201 ll) . 



3.1. Heliostat characteristics 

An heliostat is characterised by its geometry and its optical 
properties. All heliostats are assumed to be rectangular, with 
the focal equal to the slant range and made of spherical facets. 

We have a total of four parameters to describe these prop- 
erties of heliostats. All these parameters can be different for 
different groups of heliostats (what we call mixed fields): 

P01 o~h '■ Heliostat optical error. 

P02 Lh : Horizontal length of the heliostat. 

P03 L v : Vertical length of the heliostat. 

3.2. Geographic characteristics 

The location and local ground characteristics of the solar 
power plant are taken into account in the following parameters 

P05 (p '■ Latitude of the plant position. 
P06 to/v : Terrain north-south slope. 
P07 mw : Terrain east-west slope. 

3.3. Insolation and ambiente temperatures 

Direct insolation in clear days and ground level temperatures 
at day hours the 21 st of each month are part of the input data. 

Insolation ratios for each month are also input data. This ratio 
is defined as the relation between the solar energy received in a 
month and the solar energy received if all days were clear. 

4. Performance models 

We will not detail all the models used to compute the perfor- 
mance of the power plant, only to say some words and give the 
appropriate references. 

The interception efficiency is computed with an optical 
model described in (lKieral[r980albl) . 

Atmospheric losses are e stimated by the model described 
in dBiggs and Vittitoel [l976b . 

Receiver and cycle efficiency is estimated by taking into ac- 
count the total input power through the aperture, the aperture 
size, the aperture inclination and the ambient temperature. This 
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estimation is done by using a polynomial function whose coef- 
ficients have been determined by fitting data. 

Finally the total daily and yearly energies are computed by 
integrating the hourly power output. We compute total daily 
and yearly energies for clean day, and mean cloudy basis. 

5. Optimization procedure 

Within the framework described in the preceding sections we 
are in the position to start an optimization procedure. By using 
the performance models described above we can estimate the 
total yearly energy output, and with this our optimization cite- 
ria, the price of the produced energy, as a function of our eleven 
design variables ([V01-V1 1]), seven design parameters ([P01- 
P07]) and insolation data. 



E(Vi\ p a , I) 



(7) 



where v, (i = 1, . . . , 1 1) represents the optimization variables, 
p a (a — 1 , . . . , 7) are the design parameters, and / is the insola- 
tion data. 

One only needs to find the values of v, that makes Efyf, p a , I) 
maximum. Being E{vc p a ,I) a non linear function one should 
be concerned about the existence of local optima. 

To have a full control over the optimization process we have 
made some tests on typical plant configuration by using several 
state of the art algori thms. First we will use our own local op- 
timization algorithm dRamosll2008l) . designed with the particu- 
larities of solar plant optimization in mind. This is a fast algo- 
rithm that find the closest local optimum. Second we will use a 
mixture of Monte Carlo and conjugate gradient like algorithm s 



coded in the CERN MINUIT library djames and Roosl 1 1975b 



Finally we have specifically coded a genetic algorithm dRamos 



2005) commonly used to solve difficult optimization problems 
with many local optima. 

Algorithms capable of finding global optima are usually 
much slower than local algorithms. The porpoise of using 
global searchers is to find out if our preferred local algorithm 
finds the same optima as global algorithms. 

Now we will describe in detail our optimization algorithms. 

5.1. NSPOC algorithm 

Our preferred loc al algorithm is a variant of Powell's algo- 
rithm dPowelQiill). 

Starting from an approximate value of the optimal variables 
vf^ our algorithm performs a line search along the first direc- 
tion (vi) until a optimum is found with some initial precision 
h. Then the algorithm proceeds with the second direction, and 
so on, until the 1 1 variables have been explored. This process 
needs to be repeated until all the variables are fixed at their op- 
timum for a full cycle. Finally one can increase the precision of 
the optima position by repeating the process with an increased 
precision (i.e. setting h — » h/2). 

Following this algorithm we will "zig-zag" until the optima 
is found (see Fig. [SJ. One may naively think that this zig-zag 
results in a very slow convergence, but this turns out not to be 




Figure 6: "Zig-zagging" of the NSPQC algorithm to find the minimum in an 
example two dimensional function 



the case. Exploring the dimensions one by one instead of ap- 
proaching the minima along a "direct" route has the advantage 
that one does not have to fully evaluate the objective function at 
each step. For example by changing the aperture size or incli- 
nation one does not need to recompute the blocking and shad- 
owing effects, or the atmospheric efficiency. 

Moreover algorithms that approach the minima along a "di- 
rect" route needs additional information from the objective 
function in order to find the proper direction in which the op- 
tima is located. This extra information usually comes from the 
derivativ e of the objective function (i.e. the conjugate gradient 
method (IHestenes and Stiefell,ll952l) ). We do not have the pos- 
sibility of computing analytically the derivative of the objective 
function, and a numerical evaluation of the derivative severely 
worsen the convergence speed of these algorithms. There also 
exists methods that tries to avoid the '"zig-zagging" without 
using the derivative (like the original Powell method), but in 
practice this methods can easily end up searching the optima in 
a lower dimensional subspace. 

In summary, cause the re-evaluation of the objective function 
is much simpler when only one variable has changed, and be- 
cause more complicated ingredients need additional knowledge 
of the objective function to work (like the gradient) our simple 
variation of Powell's method turns out to be a fast and robust 
method for solar power plant optimization. 



5.2. MINUIT algorithm 

The MINUIT library has been widely used in high energy 
physics as well as in other fields with literally thousands of pa- 
pers based in its results. It is considered a robust minimizer. 
This is the reason we ch oose to use it here. We r efer the reader 
to the MINUIT reference (James and Roos, 1 19751) for further de- 
tails. Here we will only comment that the MINUIT library has 
several algorithms implemented. 

We first use the SEEK optimizer. This is a Monte Carlo search 
that has potentially the ability to jump over function barriers to 
find a better global optima. Th e algorithm includes a metropolis 
step (Me tropolis et al.1 Q953) which moves to a new position 
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vf m) from an old one (v;° W) ) with a probability 



id) 



(new) \ 

v; 



exp < - 



p / (new) f 



E{vf ld) ; Pa ,I) 



(8) 



We refine this optima search with the use of the MINUIT op- 
timizer MIGRAD, that is consi der MINUIT bes t local optimizer. 
It is a quasi-Newton method dDavidonl 1 1 99 lb with inexact line 
search and a stable metric updating scheme. We decided to use 
the strategy that make less use of the numerical estimates of the 
gradient of the objective function at the price of being slower. 

Regardless the potential abilities of this combination of 
MINUIT algorithms it is fair to say that about half of the MINUIT 
hackers believes that the ability of finding global minima are 
small in practical situations. 

5.3. Genetic algorithm 

To be sure that our local optimizer and MINUIT are not falling 
in local optima , we have also implemented a genetic algorithm 
dRamosll2005lPi 



Genetic algorithms (Holland, 11975b mimic the process of 
natural evolution by creating a "population" and implementing 
natural-inspired mechanisms as crossover, mutation and natural 
selection to optimize an objective function. 

The main ingredients of a genetic algorithm are 

1. A population of organisms whose fitness is given by the 
objective function. 

2. A crossover process, by which two members of the pop- 
ulation ("parents") give rise to two different members 
("childs"). 

3. A mutation process, by which organisms randomly 
change. 

4. A selection process by which the members of the popula- 
tion are chosen for later crossover and/or mutation based 
in its fitness. 

5.3.1. Population 

The members of our population have 1 1 "genes" (g,, i = 
1, . . . , 11) whose values are the 1 1 variables that code the plant 
design. In this case these variables are stored simply as a vector 
of 1 1 real numbers. 

Members of the population have a "fitness" that is given as 
the value of the objective function evaluated for its genetic con- 
tent. 



Fitness = E(gi,p a ,I) 



(9) 



We have repeated the optimization procedure with different 
population sizes in the range N, , — 30 - 100. 



4 The FORTRAN code of the genetic library is free software under the GPL 
license. 



5.3.2. Crossover and mutation 

The population is paired up, and each pair crosses with a 
probability p c (typical values in our runs are p c as 0.05). When 
a crossover occurs each of the genes of the two members of 



the population and gf produce the offspring with genetic 



,( 2 ) 



content and g. w given by 



'(2) 



„'a) 



= ngf + ii-rdgf 



Si' = Qigf + il-qdgf 



(10) 
(11) 



where r,-, q, are samples of a normal distribution A/X0.5,0.5) 
with mean 1/2 and standard deviation 1 /2. It is important here 
to remark that this choice allow interpolation as well as extrap- 
olation between parent's genet ic content, avoiding a fast "false 
convergence" (see for example lAllanach et a l. (2004)). 

After this process each gene of each member of the popula- 
tion is mutated with a probability p m (typical values in our runs 
are p m = 0.1). In this mutation process the gene value is mul- 
tiplied by (1 + r) where r ~ N(Q, 1) is a random sample of a 
standard normal distribution. 

5.3.3. Selection 

If a member is the product of a crossover or has mutated its 
fitness (value of the target function) needs to be re-computed. 

The "best" N e ute of the total population N to t organisms are 
retained from one generation to the next without change. For 
the rest we use the well known rou lette-wheel selection process 
(see for example (IGoldbergl 1 19891) and references therein) in 
which the members of the population are ordered by its fitness 
and accumulated normalized fitness values are computed (the 
accumulated fitness value is the sum of the fitness of all the 
individuals better fitted than itself). After we choose N tot —N e u te 
individuals by drawing uniform random numbers and taken the 
first individual whose accumulated normalized fitness value is 
greater than these random samples. 

The idea behind this process is that if an organism has a fit 
that is twice than other organism's fitness, it is twice as probable 
for this organism to "survive". 

This members together with the N e u te members form the 
"next generation" that is subject to the same process again. Typ- 
ically N e uie is chosen 5-10. 

5.3.4. Stop criteria 

The stop criteria is always a subtle matter in a genetic al- 
gorithm. We choose to iterate the process a fixed number of 
generations (in practice around 200), and examine if the fittest 
organism has changed in the last 20 generations. If this is the 
case we proceed to iterate for another 200 generations. 

In all situations we are sure that when we decide to stop the 
optimization process, the fitness of the elite (N e ut e ) organism 
has not significantly changed in at most 20 generations. 

We also routinely perform tests during the run to ensure that 
there exists diversity among the organisms to avoid a false fast 
convergence. The existence of diversity is controlled by com- 
puting the variance of each genes in all generations. 
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6. Impact of the variables in the plant performance 

Once one finds the value of the eleven variables that makes 
the objective function optimal, it is interesting to analyze the 
properties of this optimum. The basic idea is better understood 
in a one dimensional case (see Fig. [7). When two algorithms 
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Figure 7: One dimensional example of the minimization of a complex function. 
Algorithm 1 outputs that the minimum is at x = 1.03 with a value of the objec- 
tive function equal to 1.1808. Algorithm 2 says that the minimum is located at 
x = 1.98 with a value of the objective function equal to 1.18001. These results 
are compatible with both function types A and B, but the physical interpreta- 
tion is very different in both cases. A detailed study of the second derivative of 
the objective function will help us in discerning both cases (see text for more 
details). 

outputs the position of the minimum and the value of the objec- 
tive function in the minimum they will never give exactly the 
same answer. When the value of the objective function is very 
similar, the difference in the position of the minimum can be 
due to the fact that the minimum is very wide (function A of 
Fig. UJ, or that the function has several local minimum (func- 
tion B of Fig.|7j. These cases can be easily discerned by look- 
ing at the second derivative at the minimum. Furthermore, the 
same analysis also provides information on how the value of 
different variables affect the plant performance, helping us to 
detect "irrelevant" variables that we could have included in the 
design, and to focus in variables that are key for the plant per- 
formance. This analysis for the general case of an objective 
function of several variables will be explained in this section 



and in the Appendix A 



We will assume that the optimal value of the objective func- 
tion is achieved at v, and that the value of the objective function 
at this point is E{vj\ p a , I) = E. If we take the value of one vari- 
able away from its optimal value by an amount cr(e) the value 
of the objective function will increase. But some of the rise in 
the objective function can be compensated by tuning the other 
variables. We define the uncertainty of the variable v, (and use 
the symbol cr,(s)) with the condition that the minimum of the 
objective function with the i - th variable fixed at the value 
v; = v; + cri(e) is equal to E + s (see Fig. [5J. 

If s is chosen small, a large value of cr a (s) means that the 
variable v a plays no role in the plant performance. It does not 




o-i(£) 



Figure 8: Contour plot of the function E(Vj,p a ,I) in the Vi,V2 plane. The 
minimum of the function is located at v\ = v\ and V2 = V2 (at the center of the 
ellipse), and the value of the function in this point is E. The ellipse corresponds 
to the points that having the values of vj and vi that correspond to the point in 
the plane, have a minimum value of the objective function equal to E + e. The 
size of the ellipse defines the quantities <X[(e) and o"2(e). The inclination define 
the "correlation" between the two variables pi2- This correlation can also be 
estimated from the difference between <Xj nncr and ct(e). 



matter the value we choose for v a the objective function will 
not deviate from its optimal value. On the other hand if o- a (s) 
is small this means that any departure of v a from its optimal 
value will severely worsen the plant performance. Note that 
this uncertainty is an estimate of how wide the minimum is. 
In our one dimensional example of Fig. |7j the uncertainty of 
the variable x is about 0.7 for the function A and 0.08 for the 
function B using a value of s — 0.001. 

When missing the minimum for a variable i can be compen- 
sated by tuning the variable j we say that the variables i and j 
are correlated. They share part of the information of a plant de- 
sign. For example the tower height and the aperture inclination 
must be correlated, since we need to adjust the aperture incli- 
nation to correct an error in the tower height. We can say that 
both variables "share" the information that the aperture need to 
look directly to the field in order to increase the interception ef- 
ficiency. This correlation can be mathematically expressed by 
the correlation coefficient between the two variables py, a real 
number between -1 and 1 that measures the difference between 
the quantity cr,(e) and its would-be value if the variable j would 
have stayed fixed at its minimum value v j (the quantity cr; nner in 
Fig.©. 

A value of py close to zero means that the variables v, and 
Vj encode different information of the plant design. On the 
other hand if py is close to its extreme values py « +1, both 
variables v, and vj encode the same information, since a non 
optimal value of one of them can always be compensated by 
choosing the other appropriately. 

On one hand if a plant design variable has a very large uncer- 
tainty this variable can be dropped from the plant design, since 
it plays no role in the plant performance. On the other hand, 
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ideal plant design variables should be as uncorrelated as possi- 
ble, since otherwise the optimization algorithm will be looking 
for an optimal value of a function with an almost fiat direction. 
We are interested in computing the quantities <x,(e) and the 



correlation matrix py. As it is shown in the Appendix A when 
e is small, these quantities can be estimated by computing the 
Hessian of the objective function evaluated at the minimum 



H,j = 



d 2 E(vr, Pa j) 



(12) 



Being concise, the quantities <x,(e) are given by the diagonal 
elements of the inverse of the Hessian 



(H- l hs, 



(13) 



and the correlation coefficients are given by the normalized non 
diagonal elements of the inverse hessian 



Pi] = 



^(H- l UH- l )j 



(14) 



We emphasize here that this expressions are accurate only in the 
case that e is small enough so that all the cr,(e) are smaljf|. For 
practical purposes we can only consider the resulting cr,(e) an 
estimate of the real uncertainties, but this information is more 
than enough to obtain a qualitative understanding of how the 
different design variables impact the plant performance and to 
interpret the result of the optimization that different algorithms 
produce. 

The strategy that we follow once the optimal value of the 
design variables is found is to compute numerically the Hes- 
sian matrix, and then invert it. To compute the numerical first 
and second derivatives we use standard finite differences ex- 
pressions (we drop the non variable arguments of E for clarity): 



8E( Vi ) E( Vi + h k k) - E( Vi - h k k) 
— = hm 

dv k h^o 2h 



(15) 



d 2 E( Vi ) E( Vi + h k k) - 2E(vd + E(Vi - h k k) 

-75 = I'm (16) 

d 2 v k fc->o h 2 

d 2 E( Vi ) l f . 

a a = , lim , 7TT~ E ( Vi + hkk + h if) 
dv k dvj f><->0 4h k hj 1 



- E(vi + h k k-h jj) 

- E(Vi - h k k + h jj) 
+ E(vi - h k k + hjf)} 



(17) 



where i represents the unit vector in the direction of the variable 

Vj. 

Choosing an appropriate step size (the values for /z, in the pre- 
vious formulas) is crucial to compute derivatives numerically. 
This is a well known problem: numerical differentiation is an 
ill defined problem for finite precision arithmetic. 

We have chosen a step size hi for each direction that meets 
the following two criteria: 



5 In the almost trivial case that the objective function is quadratic the above 
mentioned expressions are exact for any value of e, but this is hardly an inter- 
esting case of study. 



1. The numerical computation of the first derivative gives a 
very small result. 

2. The numerical computation of the second derivative gives 
approximately the same result if we vary the step by a 10- 
20%. 

The first condition ensures that we are not having rounding 
errors (being at the minimum means that the first derivative 
should be zero). The second condition ensures that the com- 
putation of the second derivatives is stable. 

7. Results 

Following the steps described in the previous sections we will 
present some results about the optimal field design of some typ- 
ical solar power plants. 

For the case of a cavity receiver we will analyze the case of a 
900 heliostat power plant. This plant produces aound 15 MWe. 
We will label this plant as N900. 

For the case of the cylindrical receiver we will analyze the 
case of a 3000 heliostat field. This circular filed layout produces 
around 50 MWe. We will label this plant C3000. 

In both cases the optimization criteria consist in finding the 
cheapest price for the generated power. 

7.1. Algorithm analysis 

First we will focus in the N900 plant. In Tab. Q] we can see 
the result of the optimization procedure with each of the algo- 
rithms. As can be seen the three algorithm give a similar result 
for the value of the objective function at the minimum (quantity 
E). Nevertheless, this optimal value is achieved with different 
values of some of the variables. The value of ao is 5.4 for the 
minimum found by the NSP0C algorithm and 2.3 for the Ge- 
netic one. Since the value of the objective function is almost 
the same. This may induce to think that the objective function 
has several almost degenerate local minimums. 

But nothing could be more wrong. As our minimum analysis 
shows all these values corresponds (approximately) to one and 
the same minimum. But this minimum is very wide in some of 
the directions. This can be seen by computing the uncertainties 
needed to change the objective function value in an amount s 
for each of the directions. We will choose s = 0.001 as an 
estimate of the difference of the optimum value obtained by 
different algorithms. Clearly this small change of the objective 
function does not change the plant properties, since in the plant 
evaluation one uses approximations that lead to an estimate of 
the plant performance that is accurate with less precision than 
this value of 0.001. 

As can be seen in TabQ] the variables that show a strong 
discrepancy between the results of different algorithms corre- 
sponds to directions in the objective function in which the min- 
imum is very wide. For the case of ao the associated uncertainty 
is 2.8, meaning that the values 5.42 and 2.27 that different al- 
gorithms find correspond to the same optimum. This confirms 
that all algorithms are finding the same minimum, but that in 
this minimum the objective function has a very mild depen- 
dence with some variables. Variables that are crucial for the 



8 





NSPDC 


Genetic 


MINUIT 


cr(e = 0.001) 


E 


0.18092 


0.18113 


0.18311 


N/A 


flo 


1 A 


9 7 


j.j 


9 8 
Z.o 


ci\ x iu 


J.i J 


7 78 
J. / O 


^ 77 


U. JJ 


dg 


-y.D 


1 Q 


-J.^t 


8 1 
o. 1 


eg 


-U.UJ 


M 78 
-U. / O 


-U.Ui 


U.DO 


6 


88 
W.oo 


O 9/f 


1 8 
l.o 


1 1 
1 . 1 


6 


0.169 


0.024 


0.015 


0.097 


b x 10 3 


54.2 


32.8 


29.6 


24.0 




16.8 


19.5 


19.1 


2.5 


hj 


120. 


117. 


123. 


12. 


r 


10.78 


10.74 


10.85 


0.99 


eh 


28.6 


26.5 


38.1 


9.3 



Table 1; Comparison of the three algorithms for the N900 plant. The second, 
third and fourth column shows the results for the NSP0C, the Genetic and the 
MINUIT algorithm respectively. Finally the last column quote the values of <x(£) 
for a change in the function of e = 0.001. The first row show the value of the 
function at the minima. Subsequent rows show the value of each of the design 
variables. . The following lines shows the value of different variables at the 
minima for each algorithm. 

plant design are the ones that have a small associated uncer- 
tainty, like for example a\ or r. For this variables we can see 
that all algorithms find approximately the same value of the ob- 
jective function. 

In one sentence, the discrepancy between the value of the 
value of the variables at the minimum obtained with different 
algorithms is nicely explained once ones analyzes how strong 
is the dependence of the objective function at the minimum with 
respect to different variables. One can conclude that all the al- 
gorithms are finding the same basic plant design. 

Also we have checked that the correlations between param- 
eters are not large. The largest correlation turns out to be be- 
tween variables bo and b\, and amounts to 0.91. Note that pre- 
cisely these are the variables with a large uncertainty. Variables 
with small associated uncertainties have usually small corre- 
lations. For example in the case of r and a\ the correlation 
amounts to -0.08. 

The correlation matrix also seem to pick up important infor- 
mation of the plant design. For example the tower height hj is 
positively correlated with the aperture inclination ej (the corre- 
lation amounts to 0.1), indicating that if ones makes the tower 
higher one needs to incline more the aperture. 

In Tabf2]we can see the same information for the case of the 
C3000 plant. As we can see the conclusions are roughly the 
same. The three algorithms seem to find the same design as can 
be seen by comparing the values of the variables at the mini- 
mum found by different algorithms with their respective uncer- 
tainties. 

Comparing Tab. [1] and Tab. [2] we observe that in general the 
uncertainties are smaller for the C3000 plant. This means that 
for the design of bigger plants variables need an accurate value 
if we want to obtain a high performance. The design of solar 
power plants become more involved and difficult when the size 
of the power plant increases. This conclusion seem intuitively 
correct, since there are some variables (ej. dg that "compress" 
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0.17022 


0.16981 
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N/A 
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O 9Q 
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u. ID 


W.J J 


5 


0.109 


0.081 


0.134 


0.039 


b x 10 J 


40. 


60. 


47. 


13. 


u o 


17.8 


15.3 


15.1 


1.9 


hj 


145. 


147. 


152. 


13. 


r 


8.58 


8.66 


8.57 


0.96 


h r 


8.14 


8.15 


8.17 


0.96 



Table 2: Comparison of the three algorithms for the C3000 plant. The second, 
third and fourth column shows the results for the NSPQC, the Genetic and the 
MINUIT algorithm respectively. Finally the last column quote the values of cr(e) 
for a change in the function of £ = 0.001. The first row show the value of the 
function at the minimum (£). Subsequent rows show the value of each of the 
design variables. The following lines shows the value of different variables at 
the minimum for each algorithm. 

the field layout in the south) that seem to be crucial only for big 
circular plant designs, and more or less irrelevant for the design 
of small power plants. 



7.2. Performance analysis 

A summary of the performance of the different algorithms is 
presented in Tab. [3] 



PLANT 


Algorithm 


Calls 


Time [s] 


Time [a.u.] 




NSP0C 


1270 


17280 


1.00 


N900 


MINUIT 


3737 


86450 


5.00 




Genetic 


2381 


40539 


2.35 




NSP0C 


1872 


41460 


1.00 


C3000 


MINUIT 


4599 


102699 


2.47 




Genetic 


10398 


254829 


6.14 



Table 3: Running times of the different algorithms in our prototype plants. The 
first three columns corresponds to the 900 heliostat north field plant N900, and 
the last three columns to the 3000 heliostat cylindrical receiver plant C3000. 
The table shows the number of calls to the plant evaluation function as well 
as the time of the ran, both in seconds and relative to the running time of the 
NSP0C algorithm. As the reader can see, it took for the global algorithms (both 
MINUIT and Genetic) between 2.5 and 5 times more time to find the optimum 
plant design than the time used by NSP0C for the case of the plant N900. In the 
case of the plant C3000 the case is worse (between 4 and 7 times more time). 

The NSP0C algorithm is always faster than the global opti- 
mizers MINUIT and Genetic. In the case of the plant N900, 
MINUIT uses 5 times more time to achieve the optimum, while 
Genetic uses 2.35 more time than NSPDC. This difference is sim- 
ilar for the case of the C3000 plant. 

The two global optimizers that we have used are stochastic 
in nature, thus this running times should not be treated as ex- 
act numbers, but they are representative, and the conclusion is 
always the same: the NSPDC algorithm outperforms the global 
optimizers while giving the same results. 
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As we have said, the NSPDC code only re-compute the 
"needed" pieces of the objective function. For example block- 
ing and shadowing effects are not computed when we perform 
a line search in the direction of the aperture size. This approach 
seems to be very successful for the case of the N900 plant. On 
average each function call takes between a 25% and a 50% less 
time for the NSPDC algorithm than for the others. In the more 
involved case of the C3000 plant the NSPDC function calls still 
are faster, but with a large margin. 

8. Conclusions and perspectives 

We have presented and analyzed a method to design so- 
lar power plants. Based in the field layout done within the 
German-Spanish GAST project (1982-86) in the company IN- 
TERATOM, we have proposed a method to design solar power 
plants. This method reduces the plant design to the value of 
1 1 variables that determine the field layout, tower and receiver 
characteristics. The problem of finding optimal plant designs 
is reduced to the numerical problem of optimizing a non linear 
function. 

Due to the non linear nature of the target function and the 
high computational costs of evaluating plant performances, it 
is crucial to find a robust and fast algorithm to perform plant 
designs. Fast algorithms to optimize functions of several vari- 
ables are local in nature: they find a local optima of the target 
function, but give no information on the possible existence of 
other global optima. On the other hand global optimizers are 
much more computationally expensive and their coding is far 
more involved. 

In this work we have tackled the optimization problem with 
different algorithms, both local and global. We have shown that 
in all the cases our local NSPOC algorithm give the same re- 
sults as other more complex optimizers with significantly less 
computing effort. 

The design variables that we have choosen have enough flex- 
ibility to provide solar plant layouts with very high performance 
for a broad range of sizes: from small 10 MWe plants to big 200 
MWe plants. With our choice of variables we find that our local 
optimizer outperforms global optimizers. 

In analyzing the result of different algorithms we have also 
developed an interesting method to get information on how dif- 
ferent variables affect the performance of an optimal design. 
The method is based in analyzing the Hessian of the objective 
function at its optimum value, and gives us an estimate of how 
much the departure of a variable from its optimal value affects 
the plant performance. We have observed that circular-like big 
plants need a more accurate tuning of the variables in order to 
achieve an optimal performance. 

This information can be used to speed up the optimization 
process. It turns out to be convenient to first tune variables that 
are crucial for the optimal design, and only after worry about 
variables that have a mild impact in the plant performance. Our 
NSPOC algorithm partially profit from this information to speed 
up the convergence to the optimum. 

Moreover, this method of evaluating the impact of design 
variables in the plant performance is key to the study of the im- 



provement of plant designs. Any new proposed design should 
include a similar analysis to detect superfluous variables in the 
design and determine what are the key ingredients of the new 
design. 

This works provide the tools to address very interesting prob- 
lems, what we call parametric analysis, in which the efficiency 
of the solar power plant can be studied as a function of plant 
design parameters. Questions like Does a terrain with slope im- 
prove the plant efficiency? How much does the heliostat aspect 
ratio affect the plant efficiency? Can we build cheaper plants 
by using cheaper heliostats without losing efficiency? Some 
of these questions are currently under study by the authors, 
and prelimi nary results were pre sented in the solarPACES 20 1 1 



conference (ICrespo et all 1201 lb . 



We have also provided the basis to analyze more complex 
plant designs. An example that we consider very interest- 
ing are multi-tower layouts, in which a solar power plant is 
built with several towers, and heliostats choose which tower 
they aim based on performance (some resul ts were already pre- 
sented in the 2009 solarPACES conference (ICrespo and Ramos , 
2009)), but there are many more possibilities like multi-cavity 



receivers, or properly addressing the scalability problem in so- 
lar power plant designs: How can we build small plants that 
can later be enlarged without loosing performance? Some of 
this questions and others within the framework presented here 
are currently being studied. 

Appendix A. Computation of uncertainties and correla- 
tions. 

The techniques described in this appendix are typical in sta- 
tistical description of data. In this context the function that one 
wants to optimize is the quadratic deviation between data and 
predictions of a model (usually called^ 2 ). Thus any interested 
reader can consult any standard book on statistics for a more de- 
tailed proof of the expressions developed here. We will assume 
that we are minimizing an objective function. The expressions 
remain basically unchanged for the case of the maximization of 
a function if one changes the sign of M, (M -1 ) and the Hessian 
in the expressions that follow. 

We will start assuming that our objective function of n vari- 
ables (xj, i = 1 . . . , ri) is quadratic with a minimum located at 
Xi = Xj. The most general function with these characteristics 
can be written as 



fix) = f + 2^(* ; - Xi)Mij(Xj - xj) 



(A.l) 



where /o is the value of the function at the minimum and My is 
a symmetric positive definite matrix. To obtain cr a (s) we will 
define a new function g(Xk) (k + a) of n — 1 variables (all but the 
a - th) equal to the value of the original function with x a fixed 
at x a = x a + o-Je). 

g(Xk) = f (Xk, x a =x a + <r a (s)) 

= fo + 2u (Xi ~ x i) M iMi ~ */) 

+ 2 ^(xi - Xi)M ia o- a (s) + M aa cr 2 a (s). (A.2) 
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Taking the gradient of g(xk) and equating it to zero we obtain 
the position of the minimum value of the objective function at 
fixed x a — x a + crjs). 



dg_ 

dx k 



=> 4 = Yj ( m l )ki Mia °~ a ( e ) 



(A.3) 



where m is the minor matrix that corresponds to the element 
(a, a) (this is the matrix that results from M cutting down the 
a th column and rows), and m _I is its inverse. The value of the 
objective function at this point is given by 



8(4) = f{> 



■«(*)) 



/o + 



M aa - ^M ai (m l )..M ja 



K(s) (A.4) 



the quantity between brackets can easily be recognised as 
1 / (the inverse of the diagonal entry of the original ma- 

trix M). The condition that determines cr a (s) is that the increase 
in the function respect the value at the minimum should be s. 
So the quantity cr a (s) is given by 



o- 2 a (s) 



(AT 1 ) 

V /a 



(A.5) 



If we define cr mnel -(s) as the would-be cr a (s) without tuning 
the variable j, it is clear that a comparison between cr a (s) and 
o"inner(e) would give us information about the correlation be- 
tween variables i and j. If we define the correlation coefficient 
between variables i and j by the equation 



Pij 



(M-% 



V(M-')//(M- I ) ; 



(A.6) 



we can easily check that cr inner (e) and cr a (s) are related by 



°iUr = ( l ~ Pijrf^)- 



(A.7) 



Note that in the case p,j = the two quantities are equal. 

In the case that the objective function is not quadratic, we can 
use Taylor theorem. Close enough to the minimum any function 
is well approximated by a quadratic function. This means that 
for a general function the role of My is played by the Hessian 
evaluated at the minimum 



d 2 f(xd 



dxidxj 



+ 0(x - xf 



(A.8) 



It is worth noting that this approximation will fail if cr a (e) is 
large. In order to ensure that this does not happen one need to 
keep s sufficiently small. 
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